home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / random.i < prev    next >
Text File  |  1996-02-12  |  17KB  |  520 lines
  1. /*
  2.    RANDOM.I
  3.    Random numbers with various distributions.
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1996.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. /* Contents:
  11.  
  12.    random_x   - avoids the 2.e9 bins of the random function at a
  13.                 cost of calling random twice, can be used as a
  14.             drop-in replacement for random
  15.    random_u   - convenience routine to give uniform deviate on an
  16.                 interval other than (0,1)
  17.    random_n   - return gaussian deviate
  18.    random_ipq - arbitrary piecewise linear deviate, with optional
  19.                 power law or exponential tails
  20.    random_rej - generic implementation of rejection method, can be
  21.                 used either in conjunction with a piecewise linear
  22.         bounding function, or an arbitrary bounding function
  23.         (in the latter case, the inverse of the integral of
  24.         the bounding function must be supplied as well)
  25.  
  26.    In all cases, these routines accept a dimlist or arguments to
  27.    determine the dimensions of the returned random deviates.
  28.    Furthermore, in all cases you will get back the same sequence
  29.    of deviates no matter what the dimensionality of the calls --
  30.    for example, if at some point the call to random_n(5) returns
  31.    [.3,-.1,-.8,1.1,-.9], then if instead random_n(3) followed by
  32.    random_n(2) had been called, the return values would have been
  33.    [.3,-.1,-.8] and [1.1,-.9].
  34.  */
  35.  
  36. /* ------------------------------------------------------------------------ */
  37.  
  38. func build_dimlist(&dimlist, arg)
  39. /* DOCUMENT build_dimlist, dimlist, next_argument
  40.      build a DIMLIST, as used in the array function.  Use like this:
  41.  
  42.      func your_function(arg1, arg2, etc, dimlist, ..)
  43.      {
  44.        while (more_args()) build_dimlist, dimlist, next_arg();
  45.        ...
  46.      }
  47.  
  48.      After this, DIMLIST will be an array of the form
  49.      [#dims, dim1, dim2, ...], compounded from the multiple arguments
  50.      in the same way as the array function.  If no DIMLIST arguments
  51.      given, DIMLIST will be [] instead of [0], which will act the
  52.      same in most situations.  If that possibility is unacceptible,
  53.      you may add
  54.        if (is_void(dimlist)) dimlist= [0];
  55.      after the while loop.
  56.  */
  57. {
  58.   if (is_void(dimlist)) dimlist= [0];
  59.   else if (!dimsof(dimlist)(1)) dimlist= [1,dimlist];
  60.   if (is_void(arg)) return;
  61.   if (!dimsof(arg)(1)) {
  62.     grow, dimlist, arg;
  63.     dimlist(1)+= 1;
  64.   } else {
  65.     n= arg(1);
  66.     grow, dimlist, arg(2:1+n);
  67.     dimlist(1)+= n;
  68.   }
  69. }
  70.  
  71. /* ------------------------------------------------------------------------ */
  72.  
  73. func random_x(dimlist, ..)
  74. /* DOCUMENT random_x(dimlist)
  75.  
  76.      same as random(DIMLIST), except that random_x calls random
  77.      twice at each point, to avoid the defect that random only
  78.      can produce about 2.e9 numbers on the interval (0.,1.) (see
  79.      random for an explanation of these bins).
  80.  
  81.      You may set random=random_x to get these "better" random
  82.      numbers in every call to random.
  83.  
  84.      Unlike random, there is a chance in 1.e15 or so that random_x
  85.      may return exactly 1.0 or 0.0 (the latter may not be possible
  86.      with IEEE standard arithmetic, while the former apparently is).
  87.      Since cosmic rays are far more likely, you may as well not
  88.      worry about this.  Also, because of rounding errors, some bit
  89.      patterns may still be more likely than others, but the 0.5e-9
  90.      wide bins of random will be absent.
  91.  
  92.    SEE ALSO: random
  93.  */
  94. {
  95.   while (more_args()) build_dimlist, dimlist, next_arg();
  96.   if (is_void(dimlist)) {
  97.     dimlist= [1,2];
  98.   } else if (!dimsof(dimlist)(1)) {
  99.     dimlist= [2,2,dimlist];
  100.   } else {
  101.     dimlist= grow([dimlist(1)+1],dimlist);
  102.     dimlist(2)= 2;
  103.   }
  104.   r= random_0(dimlist);
  105.   return r(1,..) + (r(2,..)-0.5)/2147483562.
  106. }
  107.  
  108. if (is_void(random_0)) random_0= random;
  109.  
  110. /* ------------------------------------------------------------------------ */
  111.  
  112. func random_u(a, b, dimlist, ..)
  113. /* DOCUMENT random_u(a, b, dimlist)
  114.  
  115.      return uniformly distributed random numbers between A and B.
  116.      (Will never exactly equal A or B.)  The DIMLIST is as for the
  117.      array function.  Same as (b-a)*random(dimlist)+a.  If A==0,
  118.      you are better off just writing B*random(dimlist).
  119.  
  120.    SEE ALSO: random, random_x, random_n, random_ipq, random_rej
  121. */
  122. {
  123.   while (more_args()) build_dimlist, dimlist, next_arg();
  124.   return (b-a)*random(dimlist)+a;
  125. }
  126.  
  127. /* ------------------------------------------------------------------------ */
  128.  
  129. func random_n(dimlist, ..)
  130. /* DOCUMENT random_n(dimlist)
  131.  
  132.      returns an array of normally distributed random double values with
  133.      the given DIMLIST (see array function, nil for a scalar result).
  134.      The mean is 0.0 and the standard deviation is 1.0.
  135.      The algorithm follows the Box-Muller method (see Numerical Recipes
  136.      by Press et al.).
  137.  
  138.    SEE ALSO: random, random_x, random_u, random_ipq, random_rej
  139. */
  140. {
  141.   while (more_args()) build_dimlist, dimlist, next_arg();
  142.   a= array(0.0, dimlist);
  143.   scalar= !dimsof(a)(1);
  144.   if (scalar) a= [a];  /* work around long-standing Yorick bug */
  145.  
  146.   /* crucial feature of this algorithm is that the same sequence
  147.    * of random numbers is returned independent of the number requested
  148.    * each time, just like random itself */
  149.   na= numberof(a);
  150.   np= numberof(_random_n_prev);
  151.   n= (na-np+1)/2;
  152.   nr= 2*n;
  153.   if (n) {
  154.     x= random(2,n);
  155.     r= sqrt(-2.0*log(x(1,)));
  156.     theta= 2.0*pi*x(2,);
  157.     x(1,)= r*cos(theta);
  158.     x(2,)= r*sin(theta);
  159.     a(np+1:na)= x(1:na-np);
  160.   }
  161.   if (np) {
  162.     a(1)= _random_n_prev;
  163.     _random_n_prev= [];
  164.   }
  165.   if (na-np<nr) {
  166.     _random_n_prev= x(nr);
  167.   }
  168.  
  169.   if (scalar) a= a(1);  /* work around long-standing Yorick bug */
  170.   return a;
  171. }
  172.  
  173. /* storage for odd random_n values */
  174. _random_n_prev= [];
  175.  
  176. /* ------------------------------------------------------------------------ */
  177.  
  178. func random_ipq(model, dimlist, ..)
  179. /* DOCUMENT random_ipq(ipq_model, dimlist)
  180.  
  181.      returns an array of double values with the given DIMLIST (see array
  182.      function, nil for a scalar result).  The numbers are distributed
  183.      according to a piecewise linear function (possibly with power law
  184.      or exponential tails) specified by the IPQ_MODEL.  The "IPQ" stands
  185.      for "inverse piecewise quadratic", which the type of function
  186.      required to transform a uniform random deviate into the piecewise
  187.      linear distribution.  Use the ipq_setup function to compute
  188.      IPQ_MODEL.
  189.  
  190.    SEE ALSO: random, random_x, random_u, random_n, random_rej, ipq_setup
  191.  */
  192. {
  193.   while (more_args()) build_dimlist, dimlist, next_arg();
  194.   ymax= (*model(4))(1);
  195.   return ipq_compute(model, ymax*random(dimlist));
  196. }
  197.  
  198. func random_rej(target, bound, dimlist, ..)
  199. /* DOCUMENT random_rej(target_dist, ipq_model, dimlist)
  200.          or random_rej(target_dist, bounding_dist, bounding_rand, dimlist)
  201.  
  202.      returns an array of double values with the given DIMLIST (see array
  203.      function, nil for a scalar result).  The numbers are distributed
  204.      according to the TARGET_DIST function:
  205.         func target_dist(x)
  206.      returning u(x)>=0 of same number and dimensionality as x, normalized
  207.      so that the integral of target_dist(x) from -infinity to +infinity
  208.      is 1.0.  The BOUNDING_DIST function must have the same calling
  209.      sequence as TARGET_DIST:
  210.         func bounding_dist(x)
  211.      returning b(x)>=u(x) everywhere.  Since u(x) is normalized, the
  212.      integral of b(x) must be >=1.0.  Finally, BOUNDING_RAND is a
  213.      function which converts an array of uniformly distributed random
  214.      numbers on (0,1) -- as returned by random -- into an array
  215.      distributed according to BOUNDING_DIST:
  216.         func bounding_rand(uniform_x_01)
  217.      Mathematically, BOUNDING_RAND is the inverse of the integral of
  218.      BOUNDING_DIST from -infinity to x, with its input scaled to (0,1).
  219.  
  220.      If BOUNDING_DIST is not a function, then it must be an IPQ_MODEL
  221.      returned by the ipq_setup function.  In this case BOUNDING_RAND is
  222.      omitted -- ipq_compute will be used automatically.
  223.  
  224.    SEE ALSO: random, random_x, random_u, random_n, random_ipq, ipq_setup
  225.  */
  226. {
  227.   if (is_func(bound)) {
  228.     brand= dimlist;
  229.     dimlist= more_args()? next_arg() : [];
  230.   }
  231.   while (more_args()) build_dimlist, dimlist, next_arg();
  232.   if (!is_func(target) || (!is_void(brand) && !is_func(brand)) ||
  233.       (!is_func(bound) && structof(bound)!=pointer))
  234.     error, "improper calling sequence, try help,random_rej";
  235.   if (!is_func(bound)) ymax= (*bound(4))(1);
  236.  
  237.   /* build result to requested shape */
  238.   x= array(0.0, dimlist);
  239.   nreq= nx= numberof(x);
  240.   ix= 1;
  241.  
  242.   do {
  243.     /* get 25% more pairs of random numbers than nreq in order
  244.      * to allow for some to be rejected -- should actually go for
  245.      * integral(bounding_dist) times nreq, but don't know what
  246.      * that is -- could refine the estimate as each pass gets a
  247.      * better notion of the fraction rejected, but don't bother */
  248.     r= random(2, max(nreq+nreq/4,10));
  249.  
  250.     /* first get xx distributed according to bounding_rand,
  251.      * then accept according to the second random number
  252.      * continue until at least one is accepted */
  253.     for (xx=[] ; !numberof(xx) ; xx=xx(list)) {
  254.       if (is_func(bound)) {
  255.     xx= brand(r(1,..));
  256.     list= where(bound(xx)*r(2,..) <= target(xx));
  257.       } else {
  258.     xx= ipq_compute(bound, ymax*r(1,..));
  259.     list= where(ipq_function(bound,xx)*r(2,..) <= target(xx));
  260.       }
  261.     }
  262.     nxx= numberof(xx);
  263.     if (nxx>nreq) {
  264.       xx= xx(1:nreq);
  265.       nxx= nreq;
  266.     }
  267.  
  268.     nreq-= nxx;
  269.     x(ix:nx-nreq)= xx;
  270.     ix+= nxx;
  271.   } while (nreq);
  272.  
  273.   return x;
  274. }
  275.  
  276. func ipq_setup(x,u,power=,slope=)
  277. /* DOCUMENT model= ipq_setup(x, u)
  278.          or model= ipq_setup(x, u, power=[pleft,prght])
  279.      or model= ipq_setup(x, u, power=[pleft,prght], slope=[sleft,srght])
  280.  
  281.      compute a model for the ipq_compute function, which computes the
  282.      inverse of a piecewise quadratic function.  This function occurs
  283.      when computing random numbers distributed according to a piecewise
  284.      linear function.  The piecewise linear function is u(x), determined
  285.      by the discrete points X and U input to ipq_setup.  None of the
  286.      values of U may be negative, and X must be strictly increasing,
  287.      X(i)<X(i+1).
  288.  
  289.      If U(1) or U(0) is not zero, you may want to model the "tail" of
  290.      the distribution, which is not well modeled by a piecewise linear
  291.      function.  You can specify a power law tail using the power=
  292.      keyword; the left tail has initial slope (U(2)-U(1))/(X(2)-X(1)),
  293.      and decays to the left as 1/X^PLEFT.  Similarly for the right tail.
  294.      If U(1)==0, PLEFT is ignored, as is PRGHT when U(0)==0.  Use the
  295.      slope= keyword to specify an alternative value for the slope.
  296.      Note that PLEFT and PRGHT must each be greater than 1.0, and that
  297.      SLEFT>0 while SRGHT<0.  If either power is greater than or equal to
  298.      100, an exponential tail will be used.  As a convenience, you may
  299.      also specify PLEFT or PRGHT of 0 to get an exponential tail.
  300.  
  301.      Note: ipq_function(model, xp) returns the function values u(xp) at
  302.      the points xp, including the tails (if any).  ipq_compute(model, yp)
  303.      returns the xp for which (integral from -infinity to xp) of u(x)
  304.      equals yp; i.e.- the inverse of the piecewise quadratic.
  305.  
  306.    SEE ALSO: random_ipq, random_rej
  307.  */
  308. {
  309.   x= double(x);
  310.   u= double(u);
  311.   if (dimsof(x)(1)!=1 || numberof(x)<2 || dimsof(u)(1)!=1 ||
  312.       numberof(u)!=numberof(x) || anyof(u<0.)) error, "bad U or X arrays";
  313.  
  314.   /* compute the integral of u(x), starting from x(1),
  315.    * both at the given points x and at the midpoints of the intervals
  316.    * integ(u,x,xx) is the basic piecewise quadratic function */
  317.   bins= (u(zcen)*x(dif))(cum);
  318.   cens= x(pcen);  /* right shape, wrong values */
  319.   yc= integ(u,x, x(zcen));
  320.   dy= bins(dif);
  321.   /* note that these cens are constrained to lie between -1 and +1 */
  322.   cens(2:-1)= 4.*(yc-bins(1:-1))/(dy+!dy) - 2.;
  323.  
  324.   ymax= bins(0);
  325.  
  326.   if (!is_void(power)) {
  327.     if (dimsof(power)(1)!=1 || numberof(power)!=2 ||
  328.     anyof(power<=1.&power!=0.)) error, "illegal power= keyword";
  329.     if (!power(1) || !u(1)) power(1)= 100.;
  330.     if (!power(0) || !u(0)) power(0)= 100.;
  331.  
  332.     if (is_void(slope))
  333.       slope= [(u(2)-u(1))/(x(2)-x(1)), (u(0)-u(-1))/(x(0)-x(-1))];
  334.     if (dimsof(slope)(1)!=1 || numberof(slope)!=2 ||
  335.     slope(1)<0. || slope(0)>0.)
  336.       error, "illegal slope= keyword, or upward slope at endpoint";
  337.  
  338.     cens(1)= u(1)? slope(1)/u(1) : 1000./(x(2)-x(1));
  339.     cens(0)= u(0)? -slope(0)/u(1) : 1000./(x(0)-x(-1));
  340.     yi= u(1)/cens(1);
  341.     if (power(1)<100.) yi*= power(1)/(power(1)-1.);
  342.     ymax+= yi;
  343.     bins+= yi;
  344.     yi= u(0)/cens(0);
  345.     if (power(0)<100.) yi*= power(0)/(power(0)-1.);
  346.     ymax+= yi;
  347.  
  348.   } else {
  349.     power= [100.,100.];
  350.     cens(1)= 1000./(x(2)-x(1));
  351.     cens(0)= 1000./(x(0)-x(-1));
  352.   }
  353.  
  354.   parm= [ymax, power(1), power(0)];
  355.  
  356.   return [&bins, &x, &cens, &parm, &u];
  357. }
  358.  
  359. /* ------------------------------------------------------------------------ */
  360.  
  361. func ipq_compute(model, y)
  362. {
  363.   /*
  364.    * model= [&bins, &vals, &cens, &parm] where:
  365.    * bins   values of y, a piecewise quadratic function of x
  366.    * vals   values of x that go with bins
  367.    * cens   4*(yc-y0)/(y1-y0) - 2 where yc is value of y at
  368.    *             x=vals(pcen), except for first and last points
  369.    *             which are du/dx / u0 for the extrapolation model
  370.    * parm   [ymax, left_power, right_power]
  371.    *             maximum possible value of y, and
  372.    *             [left,right] powers (>1.0) of x for extrap. model
  373.    */
  374.   local bins, vals, cens, parm;
  375.   eq_nocopy, bins, *model(1);
  376.   eq_nocopy, vals, *model(2);
  377.   eq_nocopy, cens, *model(3);
  378.   eq_nocopy, parm, *model(4);
  379.  
  380.   i= digitize(y, bins);
  381.  
  382.   mask0= (i>1);
  383.   list= where(mask0);
  384.   if (numberof(list)) {
  385.     yy= y(list);
  386.     ii= i(list);
  387.     mask= (ii<=numberof(bins));
  388.     list= where(mask);
  389.     if (numberof(list)) {
  390.       /* handle piecewise quadratic part */
  391.       j= ii(list);
  392.       yb= bins(j);
  393.       xb= vals(j);
  394.       aa= cens(j);   /* 4*(yc-y0)/(y1-y0) - 2 */
  395.       j-= 1;
  396.       ya= bins(j);
  397.       xa= vals(j);
  398.       bb= 0.5*(1.+aa);
  399.       yq= (yy(list)-ya)/(yb-ya);
  400.       xq= bb+sqrt(bb*bb-aa*yq);
  401.       xq= xa + (xb-xa)*( yq/(xq+!xq) );
  402.     }
  403.  
  404.     list= where(!mask);
  405.     if (numberof(list)) {
  406.       /* handle right tail */
  407.       ymax= parm(1);
  408.       if (ymax > bins(0)) {
  409.     yy= (ymax - yy(list))/(ymax - bins(0));
  410.     xa= vals(0);
  411.     aa= cens(0);    /* du/dx / u0 */
  412.     pp= parm(2);    /* power */
  413.     yy0= yy<=0.0;
  414.     yy= max(yy,0.0)+yy0;
  415.     if (pp>=100.) xt= -log(yy)/aa;
  416.     else xt= (pp/aa) * (yy^(1./(1.-pp)) - 1.0);
  417.     xt+= xa + 1.e9*yy0*(xa-vals(1));
  418.       } else {
  419.     xt= array(vals(0), numberof(list));
  420.       }
  421.     }
  422.  
  423.     xq= merge(xq, xt, mask);
  424.   }
  425.  
  426.   list= where(!mask0);
  427.   if (numberof(list)) {
  428.     /* handle left tail */
  429.     if (bins(1)) {
  430.       yy= y(list)/bins(1);
  431.       xa= vals(1);
  432.       aa= cens(1);    /* du/dx / u0 */
  433.       pp= parm(3);    /* power */
  434.       yy0= yy<=0.0;
  435.       yy= max(yy,0.0)+yy0;
  436.       if (pp>=100.) xt= log(yy)/aa;
  437.       else xt= (pp/aa) * (1.0 - yy^(1./(1.-pp)));
  438.       xt+= xa - 1.e9*yy0*(vals(0)-xa);
  439.     } else {
  440.       xt= array(vals(1), numberof(list));
  441.     }
  442.   } else {
  443.     xt= [];
  444.   }
  445.  
  446.   return merge(xq, xt, mask0);
  447. }
  448.  
  449. func ipq_function(model, x)
  450. {
  451.   /*
  452.    * model= [&bins, &vals, &cens, &parm, &valu] where:
  453.    * bins   values of y, a piecewise quadratic function of x
  454.    * vals   values of x that go with bins
  455.    * cens   4*(yc-y0)/(y1-y0) - 2 where yc is value of y at
  456.    *             x=vals(pcen), except for first and last points
  457.    *             which are du/dx / u0 for the extrapolation model
  458.    * parm   [ymax, left_power, right_power]
  459.    *             maximum possible value of y, and
  460.    *             [left,right] powers (>1.0) of x for extrap. model
  461.    * valu   values of u that go with x
  462.    */
  463.   local vals, cens, parm, valu;
  464.   eq_nocopy, vals, *model(2);
  465.   eq_nocopy, cens, *model(3);
  466.   eq_nocopy, parm, *model(4);
  467.   eq_nocopy, valu, *model(5);
  468.  
  469.   i= digitize(x, vals);
  470.  
  471.   mask0= (i>1);
  472.   list= where(mask0);
  473.   if (numberof(list)) {
  474.     xx= x(list);
  475.     ii= i(list);
  476.     mask= (ii<=numberof(vals));
  477.     list= where(mask);
  478.     if (numberof(list)) {
  479.       /* handle piecewise linear part */
  480.       uq= interp(valu, vals, xx(list));
  481.     }
  482.  
  483.     list= where(!mask);
  484.     if (numberof(list)) {
  485.       /* handle right tail */
  486.       if (valu(0)) {
  487.     xx= xx(list) - vals(0);
  488.     aa= cens(0);    /* du/dx / u0 */
  489.     pp= parm(2);    /* power */
  490.     if (pp>=100.) ut= valu(0)*exp(-aa*xx);
  491.     else ut= valu(0) / (1. + (aa/pp)*xx)^pp;
  492.       } else {
  493.     ut= array(0.0, numberof(list));
  494.       }
  495.     }
  496.  
  497.     uq= merge(uq, ut, mask);
  498.   }
  499.  
  500.   list= where(!mask0);
  501.   if (numberof(list)) {
  502.     /* handle left tail */
  503.     if (valu(1)) {
  504.       xx= vals(1) - x(list);
  505.       aa= cens(1);    /* du/dx / u0 */
  506.       pp= parm(3);    /* power */
  507.       if (pp>=100.) ut= valu(1)*exp(-aa*xx);
  508.       else ut= valu(1) / (1. + (aa/pp)*xx)^pp;
  509.     } else {
  510.       ut= array(0.0, numberof(list));
  511.     }
  512.   } else {
  513.     ut= [];
  514.   }
  515.  
  516.   return merge(uq, ut, mask0);
  517. }
  518.  
  519. /* ------------------------------------------------------------------------ */
  520.